R Geospatial Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Execute chunks of code by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(tidyr)
library(ggplot2)

library(rgdal)

Loading required package: sp
rgdal: version: 1.4-6, (SVN revision 841)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28
 Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/gdal
 GDAL binary built with GEOS: FALSE 
 Loaded PROJ.4 runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
 Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/proj
 Linking to sp version: 1.3-1 

library(sp)
library(raster)

Attaching package: ‘raster’

The following object is masked from ‘package:tidyr’:

    extract

The following object is masked from ‘package:dplyr’:

    select

library(sf)

Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0

list.files()

[1] "data"                          "r_geospatial_notebook.nb.html"
[3] "r_geospatial_notebook.Rmd"     "README.md"                    

GDALinfo("./data/HARV_dsmCrop.tif")

rows        1367 
columns     1697 
bands       1 
lower left origin.x        731453 
lower left origin.y        4712471 
res.x       1 
res.y       1 
ysign       -1 
oblique.x   0 
oblique.y   0 
driver      GTiff 
projection  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
file        ./data/HARV_dsmCrop.tif 
apparent band summary:
   GDType hasNoDataValue NoDataValue blockSize1 blockSize2
1 Float64           TRUE       -9999          1       1697
apparent band statistics:
    Bmin   Bmax    Bmean      Bsd
1 305.07 416.07 359.8531 17.83169
Metadata:
AREA_OR_POINT=Area 

dsm <- raster("./data/HARV_dsmCrop.tif")
dsm

class       : RasterLayer 
dimensions  : 1367, 1697, 2319799  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 731453, 733150, 4712471, 4713838  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : /Users/francesdavenport/Desktop/r-geospatial-100720-master/data/HARV_dsmCrop.tif 
names       : HARV_dsmCrop 
values      : 305.07, 416.07  (min, max)

class(dsm)

[1] "RasterLayer"
attr(,"package")
[1] "raster"

plot(dsm)

summary(dsm)

summary is an estimate based on a sample of 1e+05 cells (4.31% of all cells)

        HARV_dsmCrop
Min.          305.56
1st Qu.       345.67
Median        359.70
3rd Qu.       374.24
Max.          415.40
NA's            0.00

cellStats(dsm, stat='mean')

[1] 359.8531

dsm[1:8, 1:8]

 [1] 408.77 408.23 406.53 406.61 409.06 409.45 413.97 415.77 407.05 406.62 404.98
[12] 405.01 407.44 407.17 413.21 414.12 407.06 406.03 403.55 404.30 408.55 411.42
[23] 409.59 409.82 405.70 406.32 404.16 403.81 404.19 403.71 405.83 408.20 404.57
[34] 412.16 396.06 392.26 401.96 402.64 397.24 396.84 407.40 413.08 412.37 407.52
[45] 407.39 405.97 399.11 403.42 411.96 413.86 413.08 414.37 410.71 406.57 399.31
[56] 402.24 412.48 414.86 414.68 415.77 414.16 414.14 401.23 413.67

crs(dsm)

CRS arguments:
 +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
+towgs84=0,0,0 

dsm_df <- as.data.frame(dsm, xy=TRUE)
head(dsm_df)

         x       y HARV_dsmCrop
1 731453.5 4713838       408.77
2 731454.5 4713838       408.23
3 731455.5 4713838       406.53
4 731456.5 4713838       406.61
5 731457.5 4713838       409.06
6 731458.5 4713838       409.45

nrow(dsm)

[1] 1367

ncol(dsm)

[1] 1697

nrow(dsm_df)

[1] 2319799

1367*1697

[1] 2319799

ggplot() + 
  geom_raster(data = dsm_df, aes(x=x, y=y, fill=HARV_dsmCrop)) + 
  scale_fill_viridis_c()

ggplot() + 
  geom_histogram(data = dsm_df, aes(HARV_dsmCrop), bins=40)

dsm_df <- dsm_df %>%
  mutate(elevation_bin = cut(HARV_dsmCrop, breaks = c(300, 340, 380, 420)))

ggplot() + 
  geom_bar(data = dsm_df, aes(elevation_bin))

dtm <- raster("./data/HARV_dtmCrop.tif")
plot(dtm)

canopy_height <- dsm - dtm
plot(canopy_height)

mean_dtm <- cellStats(dtm, stat='mean')
mean_dtm

[1] 344.8979

dtm_anom <- dtm - mean_dtm
plot(dtm_anom)

writeRaster(canopy_height, "./data/canopy_height.tif", overwrite=TRUE)

harv_rgb <- stack("./data/HARV_RGB_Ortho.tif")
harv_rgb

class       : RasterStack 
dimensions  : 2317, 3073, 7120141, 3  (nrow, ncol, ncell, nlayers)
resolution  : 0.25, 0.25  (x, y)
extent      : 731998.5, 732766.8, 4712956, 4713536  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
names       : HARV_RGB_Ortho.1, HARV_RGB_Ortho.2, HARV_RGB_Ortho.3 
min values  :                0,                0,                0 
max values  :              255,              255,              255 

plot(harv_rgb)

plotRGB(harv_rgb, r = 1, g = 2, b = 3)

aoi_sp <- shapefile("./data/HarClip_UTMZ18.shp")
aoi_sp

class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : 732128, 732251.1, 4713209, 4713359  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 1
names       : id 
value       :  1 

plot(aoi_sp)

aoi_sp[1,]

class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : 732128, 732251.1, 4713209, 4713359  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 1
names       : id 
value       :  1 

aoi <- st_read("./data/HarClip_UTMZ18.shp")

Reading layer `HarClip_UTMZ18' from data source `/Users/francesdavenport/Desktop/r-geospatial-100720-master/data/HarClip_UTMZ18.shp' using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 732128 ymin: 4713209 xmax: 732251.1 ymax: 4713359
epsg (SRID):    32618
proj4string:    +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs

plot(aoi)

head(aoi)

Simple feature collection with 1 feature and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 732128 ymin: 4713209 xmax: 732251.1 ymax: 4713359
epsg (SRID):    32618
proj4string:    +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
  id                       geometry
1  1 POLYGON ((732128 4713359, 7...

ggplot() +
  geom_sf(data = aoi)

ne_state <- st_read("./data/Boundary-US-State-NEast.shp")

Reading layer `Boundary-US-State-NEast' from data source `/Users/francesdavenport/Desktop/r-geospatial-100720-master/data/Boundary-US-State-NEast.shp' using driver `ESRI Shapefile'
Simple feature collection with 12 features and 9 fields
geometry type:  MULTIPOLYGON
dimension:      XYZ
bbox:           xmin: -80.51989 ymin: 37.91685 xmax: -66.9499 ymax: 47.45716
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

plot(ne_state)

names(ne_state)

 [1] "STATEFP"  "STATENS"  "AFFGEOID" "GEOID"    "STUSPS"   "NAME"     "LSAD"    
 [8] "ALAND"    "AWATER"   "geometry"

head(ne_state)

Simple feature collection with 6 features and 9 fields
geometry type:  MULTIPOLYGON
dimension:      XYZ
bbox:           xmin: -79.76212 ymin: 37.91685 xmax: -66.9499 ymax: 47.45716
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
  STATEFP  STATENS    AFFGEOID GEOID STUSPS                 NAME LSAD        ALAND
1      11 01702382 0400000US11    11     DC District of Columbia   00    158350578
2      24 01714934 0400000US24    24     MD             Maryland   00  25147575220
3      36 01779796 0400000US36    36     NY             New York   00 122054577774
4      09 01779780 0400000US09    09     CT          Connecticut   00  12542396439
5      23 01779787 0400000US23    23     ME                Maine   00  79885244456
6      25 00606926 0400000US25    25     MA        Massachusetts   00  20203936287
       AWATER                       geometry
1    18633500 MULTIPOLYGON Z (((-77.11976...
2  6983455225 MULTIPOLYGON Z (((-76.04621...
3 19242052501 MULTIPOLYGON Z (((-72.01893...
4  1814978794 MULTIPOLYGON Z (((-73.69594...
5 11749091526 MULTIPOLYGON Z (((-68.92401...
6  7131805598 MULTIPOLYGON Z (((-70.27553...

ggplot() + 
  geom_sf(data = ne_state, aes(fill = NAME))

ggplot() +
  geom_sf(data = ne_state, aes(fill = AWATER))

mass <- filter(ne_state, NAME == "Massachusetts")
ggplot() + 
  geom_sf(data = mass)

mass

Simple feature collection with 1 feature and 9 fields
geometry type:  MULTIPOLYGON
dimension:      XYZ
bbox:           xmin: -73.50814 ymin: 41.23796 xmax: -69.92826 ymax: 42.88459
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
  STATEFP  STATENS    AFFGEOID GEOID STUSPS          NAME LSAD       ALAND     AWATER
1      25 00606926 0400000US25    25     MA Massachusetts   00 20203936287 7131805598
                        geometry
1 MULTIPOLYGON Z (((-70.27553...

roads <- st_read("./data/HARV_roads.shp")

Reading layer `HARV_roads' from data source `/Users/francesdavenport/Desktop/r-geospatial-100720-master/data/HARV_roads.shp' using driver `ESRI Shapefile'
Simple feature collection with 13 features and 15 fields
geometry type:  MULTILINESTRING
dimension:      XY
bbox:           xmin: 730741.2 ymin: 4711942 xmax: 733295.5 ymax: 4714260
epsg (SRID):    32618
proj4string:    +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs

names(roads)

 [1] "OBJECTID_1" "OBJECTID"   "TYPE"       "NOTES"      "MISCNOTES"  "RULEID"    
 [7] "MAPLABEL"   "SHAPE_LENG" "LABEL"      "BIKEHORSE"  "RESVEHICLE" "RECMAP"    
[13] "Shape_Le_1" "ResVehic_1" "BicyclesHo" "geometry"  

head(roads)

Simple feature collection with 6 features and 15 fields
geometry type:  MULTILINESTRING
dimension:      XY
bbox:           xmin: 730741.2 ymin: 4712685 xmax: 732232.3 ymax: 4713726
epsg (SRID):    32618
proj4string:    +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
  OBJECTID_1 OBJECTID       TYPE             NOTES MISCNOTES RULEID          MAPLABEL
1         14       48 woods road Locust Opening Rd      <NA>      5 Locust Opening Rd
2         40       91   footpath              <NA>      <NA>      6              <NA>
3         41      106   footpath              <NA>      <NA>      6              <NA>
4        211      279 stone wall              <NA>      <NA>      1              <NA>
5        212      280 stone wall              <NA>      <NA>      1              <NA>
6        213      281 stone wall              <NA>      <NA>      1              <NA>
  SHAPE_LENG             LABEL BIKEHORSE RESVEHICLE RECMAP Shape_Le_1
1 1297.35706 Locust Opening Rd         Y         R1      Y 1297.10617
2  146.29984              <NA>         Y         R1      Y  146.29983
3  676.71804              <NA>         Y         R2      Y  676.71807
4  231.78957              <NA>      <NA>       <NA>   <NA>  231.78962
5   45.50864              <NA>      <NA>       <NA>   <NA>   45.50859
6  198.39043              <NA>      <NA>       <NA>   <NA>  198.39041
                             ResVehic_1                  BicyclesHo
1    R1 - All Research Vehicles Allowed Bicycles and Horses Allowed
2    R1 - All Research Vehicles Allowed Bicycles and Horses Allowed
3 R2 - 4WD/High Clearance Vehicles Only Bicycles and Horses Allowed
4                                  <NA>                        <NA>
5                                  <NA>                        <NA>
6                                  <NA>                        <NA>
                        geometry
1 MULTILINESTRING ((730819.2 ...
2 MULTILINESTRING ((732040.2 ...
3 MULTILINESTRING ((732057 47...
4 MULTILINESTRING ((731903.6 ...
5 MULTILINESTRING ((732039.1 ...
6 MULTILINESTRING ((732056.2 ...

roads[1,]

Simple feature collection with 1 feature and 15 fields
geometry type:  MULTILINESTRING
dimension:      XY
bbox:           xmin: 730741.2 ymin: 4712685 xmax: 731960.4 ymax: 4713205
epsg (SRID):    32618
proj4string:    +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
  OBJECTID_1 OBJECTID       TYPE             NOTES MISCNOTES RULEID          MAPLABEL
1         14       48 woods road Locust Opening Rd      <NA>      5 Locust Opening Rd
  SHAPE_LENG             LABEL BIKEHORSE RESVEHICLE RECMAP Shape_Le_1
1   1297.357 Locust Opening Rd         Y         R1      Y   1297.106
                          ResVehic_1                  BicyclesHo
1 R1 - All Research Vehicles Allowed Bicycles and Horses Allowed
                        geometry
1 MULTILINESTRING ((730819.2 ...

roads$TYPE

 [1] woods road footpath   footpath   stone wall stone wall stone wall stone wall
 [8] stone wall stone wall boardwalk  woods road woods road woods road
Levels: boardwalk footpath stone wall woods road

roads[,"LABEL"]

Simple feature collection with 13 features and 1 field
geometry type:  MULTILINESTRING
dimension:      XY
bbox:           xmin: 730741.2 ymin: 4711942 xmax: 733295.5 ymax: 4714260
epsg (SRID):    32618
proj4string:    +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
First 10 features:
               LABEL                       geometry
1  Locust Opening Rd MULTILINESTRING ((730819.2 ...
2               <NA> MULTILINESTRING ((732040.2 ...
3               <NA> MULTILINESTRING ((732057 47...
4               <NA> MULTILINESTRING ((731903.6 ...
5               <NA> MULTILINESTRING ((732039.1 ...
6               <NA> MULTILINESTRING ((732056.2 ...
7               <NA> MULTILINESTRING ((731964 47...
8               <NA> MULTILINESTRING ((732105.2 ...
9               <NA> MULTILINESTRING ((732222.9 ...
10              <NA> MULTILINESTRING ((732153.8 ...

st_crs(roads)

Coordinate Reference System:
  EPSG: 32618 
  proj4string: "+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs"

st_geometry(roads)

Geometry set for 13 features 
geometry type:  MULTILINESTRING
dimension:      XY
bbox:           xmin: 730741.2 ymin: 4711942 xmax: 733295.5 ymax: 4714260
epsg (SRID):    32618
proj4string:    +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
First 5 geometries:

MULTILINESTRING ((730819.2 4712714, 730812.1 47...
MULTILINESTRING ((732040.2 4713241, 732033.7 47...
MULTILINESTRING ((732057 4713726, 732053.8 4713...
MULTILINESTRING ((731903.6 4713197, 731882.6 47...
MULTILINESTRING ((732039.1 4713299, 732051.9 47...

st_coordinates(roads)

              X       Y L1 L2
  [1,] 730819.2 4712714  1  1
  [2,] 730812.1 4712714  1  1
  [3,] 730809.8 4712714  1  1
  [4,] 730808.7 4712713  1  1
  [5,] 730807.1 4712712  1  1
  [6,] 730805.4 4712710  1  1
  [7,] 730798.7 4712700  1  1
  [8,] 730797.2 4712699  1  1
  [9,] 730795.7 4712697  1  1
 [10,] 730793.1 4712695  1  1
 [11,] 730790.4 4712693  1  1
 [12,] 730787.6 4712692  1  1
 [13,] 730784.8 4712690  1  1
 [14,] 730781.7 4712689  1  1
 [15,] 730778.6 4712688  1  1
 [16,] 730775.4 4712687  1  1
 [17,] 730772.0 4712686  1  1
 [18,] 730768.6 4712686  1  1
 [19,] 730765.0 4712685  1  1
 [20,] 730757.5 4712685  1  1
 [21,] 730749.6 4712686  1  1
 [22,] 730741.2 4712687  1  1
 [23,] 731034.8 4712792  2  1
 [24,] 731042.7 4712794  2  1
 [25,] 731050.4 4712795  2  1
 [26,] 731058.0 4712798  2  1
 [27,] 731065.4 4712800  2  1
 [28,] 731075.9 4712804  2  1
 [29,] 731083.0 4712807  2  1
 [30,] 731089.2 4712810  2  1
 [31,] 731094.4 4712814  2  1
 [32,] 731097.9 4712817  2  1
 [33,] 731101.6 4712821  2  1
 [34,] 731105.4 4712825  2  1
 [35,] 731111.4 4712833  2  1
 [36,] 731115.7 4712838  2  1
 [37,] 731120.2 4712845  2  1
 [38,] 731133.2 4712865  2  1
 [39,] 731135.8 4712870  2  1
 [40,] 731143.3 4712883  2  1
 [41,] 731146.0 4712888  2  1
 [42,] 731148.6 4712891  2  1
 [43,] 731152.1 4712895  2  1
 [44,] 731154.9 4712898  2  1
 [45,] 731157.0 4712900  2  1
 [46,] 731159.4 4712902  2  1
 [47,] 731162.2 4712903  2  1
 [48,] 731165.2 4712904  2  1
 [49,] 731171.6 4712906  2  1
 [50,] 731184.8 4712909  2  1
 [51,] 731190.2 4712910  2  1
 [52,] 731195.5 4712912  2  1
 [53,] 731199.9 4712914  2  1
 [54,] 731204.3 4712917  2  1
 [55,] 731208.9 4712920  2  1
 [56,] 731212.8 4712924  2  1
 [57,] 731216.7 4712928  2  1
 [58,] 731220.8 4712932  2  1
 [59,] 731232.3 4712946  2  1
 [60,] 731236.2 4712950  2  1
 [61,] 731240.0 4712954  2  1
 [62,] 731247.2 4712960  2  1
 [63,] 731249.7 4712962  2  1
 [64,] 731252.0 4712963  2  1
 [65,] 731255.0 4712964  2  1
 [66,] 731258.3 4712966  2  1
 [67,] 731261.9 4712967  2  1
 [68,] 731265.8 4712968  2  1
 [69,] 731269.1 4712968  2  1
 [70,] 731272.9 4712968  2  1
 [71,] 731285.3 4712967  2  1
 [72,] 731289.8 4712967  2  1
 [73,] 731298.3 4712968  2  1
 [74,] 731305.6 4712968  2  1
 [75,] 731312.2 4712969  2  1
 [76,] 731318.3 4712970  2  1
 [77,] 731325.9 4712971  2  1
 [78,] 731332.4 4712973  2  1
 [79,] 731338.2 4712974  2  1
 [80,] 731343.5 4712976  2  1
 [81,] 731350.0 4712979  2  1
 [82,] 731355.6 4712981  2  1
 [83,] 731360.6 4712984  2  1
 [84,] 731365.0 4712987  2  1
 [85,] 731367.5 4712989  2  1
 [86,] 731370.3 4712991  2  1
 [87,] 731378.8 4713000  2  1
 [88,] 731382.8 4713004  2  1
 [89,] 731393.1 4713016  2  1
 [90,] 731396.5 4713019  2  1
 [91,] 731399.8 4713022  2  1
 [92,] 731406.3 4713026  2  1
 [93,] 731414.7 4713030  2  1
 [94,] 731417.9 4713032  2  1
 [95,] 731421.6 4713033  2  1
 [96,] 731434.2 4713036  2  1
 [97,] 731438.6 4713038  2  1
 [98,] 731442.7 4713039  2  1
 [99,] 731446.1 4713041  2  1
[100,] 731453.4 4713046  2  1
[101,] 731467.1 4713056  2  1
[102,] 731476.4 4713062  2  1
[103,] 731485.1 4713067  2  1
[104,] 731493.3 4713070  2  1
[105,] 731496.8 4713072  2  1
[106,] 731500.8 4713073  2  1
[107,] 731514.3 4713076  2  1
[108,] 731519.1 4713077  2  1
[109,] 731523.6 4713078  2  1
[110,] 731527.5 4713080  2  1
[111,] 731532.7 4713083  2  1
[112,] 731538.2 4713087  2  1
[113,] 731547.9 4713095  2  1
[114,] 731552.9 4713098  2  1
[115,] 731559.9 4713102  2  1
[116,] 731567.8 4713106  2  1
[117,] 731577.0 4713110  2  1
[118,] 731586.9 4713113  2  1
[119,] 731604.9 4713118  2  1
[120,] 731614.1 4713121  2  1
[121,] 731619.3 4713122  2  1
[122,] 731623.3 4713123  2  1
[123,] 731633.7 4713123  2  1
[124,] 731644.4 4713124  2  1
[125,] 731652.4 4713125  2  1
[126,] 731659.5 4713127  2  1
[127,] 731665.8 4713129  2  1
[128,] 731674.2 4713132  2  1
[129,] 731677.4 4713133  2  1
[130,] 731680.2 4713135  2  1
[131,] 731682.4 4713137  2  1
[132,] 731687.1 4713141  2  1
[133,] 731689.5 4713142  2  1
[134,] 731692.5 4713144  2  1
[135,] 731696.4 4713145  2  1
[136,] 731706.3 4713148  2  1
[137,] 731711.4 4713150  2  1
[138,] 731719.4 4713154  2  1
[139,] 731725.6 4713157  2  1
[140,] 731730.1 4713159  2  1
[141,] 731734.6 4713162  2  1
[142,] 731746.1 4713171  2  1
[143,] 731758.1 4713178  2  1
[144,] 731768.9 4713185  2  1
[145,] 731779.2 4713190  2  1
[146,] 731789.1 4713195  2  1
[147,] 731796.6 4713197  2  1
[148,] 731803.6 4713200  2  1
[149,] 731810.9 4713202  2  1
[150,] 731815.8 4713203  2  1
[151,] 731820.9 4713204  2  1
[152,] 731827.8 4713205  2  1
[153,] 731834.5 4713205  2  1
[154,] 731841.9 4713204  2  1
[155,] 731848.7 4713204  2  1
[156,] 731852.9 4713204  2  1
[157,] 731857.3 4713202  2  1
[158,] 731865.1 4713199  2  1
[159,] 731880.8 4713194  2  1
[160,] 731891.1 4713190  2  1
[161,] 731900.3 4713186  2  1
[162,] 731908.3 4713182  2  1
[163,] 731915.0 4713178  2  1
[164,] 731920.8 4713174  2  1
[165,] 731926.0 4713170  2  1
[166,] 731930.7 4713166  2  1
[167,] 731936.2 4713162  2  1
[168,] 731939.8 4713158  2  1
[169,] 731942.9 4713155  2  1
[170,] 731945.3 4713152  2  1
[171,] 731947.9 4713148  2  1
[172,] 731950.2 4713143  2  1
[173,] 731952.3 4713138  2  1
[174,] 731954.5 4713131  2  1
[175,] 731955.7 4713127  2  1
[176,] 731956.7 4713122  2  1
[177,] 731958.7 4713105  2  1
[178,] 731959.9 4713094  2  1
[179,] 731960.4 4713086  2  1
[180,] 731960.3 4713079  2  1
[181,] 731959.7 4713069  2  1
[182,] 731959.2 4713062  2  1
[183,] 731957.4 4713047  2  1
[184,] 731956.9 4713040  2  1
[185,] 731956.5 4713030  2  1
[186,] 731956.4 4713025  2  1
[187,] 731956.8 4713019  2  1
[188,] 731957.5 4713015  2  1
[189,] 731958.6 4713010  2  1
[190,] 731960.1 4713006  2  1
[191,] 732040.2 4713241  1  2
[192,] 732033.7 4713235  1  2
[193,] 732022.1 4713230  1  2
[194,] 732013.8 4713222  1  2
[195,] 732009.0 4713212  1  2
[196,] 732004.9 4713200  1  2
[197,] 732002.9 4713187  1  2
[198,] 731999.3 4713175  1  2
[199,] 731993.5 4713164  1  2
[200,] 731983.4 4713149  1  2
[201,] 731971.3 4713141  1  2
[202,] 731954.5 4713131  1  2
[203,] 732057.0 4713726  1  3
[204,] 732053.8 4713719  1  3
[205,] 732053.9 4713711  1  3
[206,] 732054.9 4713703  1  3
[207,] 732058.7 4713694  1  3
[208,] 732064.4 4713680  1  3
[209,] 732071.0 4713673  1  3
[210,] 732076.7 4713659  1  3
[211,] 732078.8 4713644  1  3
[212,] 732079.9 4713635  1  3
[213,] 732081.8 4713626  1  3
[214,] 732102.2 4713617  1  3
[215,] 732143.9 4713602  1  3
[216,] 732167.2 4713590  1  3
[217,] 732200.2 4713560  1  3
[218,] 732221.3 4713531  1  3
[219,] 732225.2 4713499  1  3
[220,] 732219.1 4713477  1  3
[221,] 732213.8 4713460  1  3
[222,] 732211.2 4713446  1  3
[223,] 732226.9 4713416  1  3
[224,] 732232.3 4713401  1  3
[225,] 732229.3 4713387  1  3
[226,] 732206.2 4713354  1  3
[227,] 732157.0 4713309  1  3
[228,] 732138.3 4713290  1  3
[229,] 732124.1 4713282  1  3
[230,] 732108.8 4713278  1  3
[231,] 732095.5 4713275  1  3
[232,] 732079.2 4713259  1  3
[233,] 732061.1 4713245  1  3
[234,] 732046.4 4713242  1  3
[235,] 732040.2 4713241  1  3
[236,] 731903.6 4713197  1  4
[237,] 731882.6 4713287  1  4
[238,] 732021.8 4713296  1  4
[239,] 732039.1 4713299  1  5
[240,] 732051.9 4713255  1  5
[241,] 732056.2 4713235  1  6
[242,] 732105.2 4713043  1  6
[243,] 731964.0 4713019  1  7
[244,] 732105.2 4713043  1  7
[245,] 732105.2 4713043  1  8
[246,] 732194.2 4713059  1  8
[247,] 732222.9 4713063  1  9
[248,] 732258.2 4713070  1  9
[249,] 732153.8 4713305  1 10
[250,] 732159.8 4713303  1 10
 [ reached getOption("max.print") -- omitted 109 rows ]

tower <- st_read("./data/HARVtower_UTM18N.shp")

Reading layer `HARVtower_UTM18N' from data source `/Users/francesdavenport/Desktop/r-geospatial-100720-master/data/HARVtower_UTM18N.shp' using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 14 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 732183.2 ymin: 4713265 xmax: 732183.2 ymax: 4713265
epsg (SRID):    32618
proj4string:    +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs

ggplot() + 
  geom_sf(data = aoi, fill = "green", alpha = 0.5) + 
  geom_sf(data = roads) + 
  geom_sf(data = tower, col = "red", size = 10, shape = "*")

ggplot() + 
  geom_sf(data = roads, aes(col = LABEL, linetype = TYPE))

tower_buffer <- st_buffer(tower, 200)

ggplot()+ 
  geom_sf(data = tower) + 
  geom_sf(data = tower_buffer, fill = "blue", alpha = 0.3)

crs(canopy_height)

CRS arguments:
 +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 

st_crs(tower_buffer)

Coordinate Reference System:
  EPSG: 32618 
  proj4string: "+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs"

canopy_masked <- mask(canopy_height, tower_buffer, inverse=TRUE)
plot(canopy_masked)

st_bbox(tower_buffer)

     xmin      ymin      xmax      ymax 
 731983.2 4713065.0  732383.2 4713465.0 

canopy_masked <- crop(canopy_masked, extent(731900, 732500, 4713000, 4713500))
plot(canopy_masked)

cellStats(canopy_masked, stat = "mean")

[1] 19.36194

canopy_extracted_values <- raster::extract(canopy_height, tower_buffer)

canopy_extracted_values <- unlist(canopy_extracted_values)
hist(canopy_extracted_values)

ndvi_files <- list.files("./data/NDVI", full.names=TRUE, pattern = ".tif$")
ndvi_files

 [1] "./data/NDVI/005_HARV_ndvi_crop.tif" "./data/NDVI/037_HARV_ndvi_crop.tif"
 [3] "./data/NDVI/085_HARV_ndvi_crop.tif" "./data/NDVI/133_HARV_ndvi_crop.tif"
 [5] "./data/NDVI/181_HARV_ndvi_crop.tif" "./data/NDVI/197_HARV_ndvi_crop.tif"
 [7] "./data/NDVI/213_HARV_ndvi_crop.tif" "./data/NDVI/229_HARV_ndvi_crop.tif"
 [9] "./data/NDVI/245_HARV_ndvi_crop.tif" "./data/NDVI/261_HARV_ndvi_crop.tif"
[11] "./data/NDVI/277_HARV_ndvi_crop.tif" "./data/NDVI/293_HARV_ndvi_crop.tif"
[13] "./data/NDVI/309_HARV_ndvi_crop.tif"

ndvi <- stack(ndvi_files)
ndvi

class       : RasterStack 
dimensions  : 5, 4, 20, 13  (nrow, ncol, ncell, nlayers)
resolution  : 30, 30  (x, y)
extent      : 239415, 239535, 4714215, 4714365  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=19 +ellps=WGS84 +units=m +no_defs 
names       : X005_HARV_ndvi_crop, X037_HARV_ndvi_crop, X085_HARV_ndvi_crop, X133_HARV_ndvi_crop, X181_HARV_ndvi_crop, X197_HARV_ndvi_crop, X213_HARV_ndvi_crop, X229_HARV_ndvi_crop, X245_HARV_ndvi_crop, X261_HARV_ndvi_crop, X277_HARV_ndvi_crop, X293_HARV_ndvi_crop, X309_HARV_ndvi_crop 
min values  :                2732,                1534,                1917,                5669,                8685,                8768,                8633,                8686,                8297,                7491,                 277,                 484,                4396 
max values  :                5545,                3736,                3600,                6394,                8903,                9140,                8945,                8967,                8651,                8607,                 366,                 661,                6359 

ndvi_Z18 <- projectRaster(ndvi, crs = crs(dtm))

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

ndvi_Z18

class       : RasterBrick 
dimensions  : 5, 4, 20, 13  (nrow, ncol, ncell, nlayers)
resolution  : 29.9, 29.9  (x, y)
extent      : 732135.1, 732254.7, 4713216, 4713365  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : X005_HARV_ndvi_crop, X037_HARV_ndvi_crop, X085_HARV_ndvi_crop, X133_HARV_ndvi_crop, X181_HARV_ndvi_crop, X197_HARV_ndvi_crop, X213_HARV_ndvi_crop, X229_HARV_ndvi_crop, X245_HARV_ndvi_crop, X261_HARV_ndvi_crop, X277_HARV_ndvi_crop, X293_HARV_ndvi_crop, X309_HARV_ndvi_crop 
min values  :           2746.0677,           1649.4942,           2037.4145,           5752.6465,           8673.7318,           8813.5836,           8583.8179,           8711.7698,           8278.1903,           7542.8330,            295.9601,            461.3890,           4732.6922 
max values  :           5715.4455,           3872.7388,           3724.4348,           6445.5631,           8869.5036,           9138.6236,           8963.2645,           8944.5587,           8610.9731,           8405.7861,            368.3793,            630.1992,           6470.4486 

plot(ndvi_Z18, 1)

cellStats(ndvi_Z18, stat="mean")

X005_HARV_ndvi_crop X037_HARV_ndvi_crop X085_HARV_ndvi_crop X133_HARV_ndvi_crop 
          3709.2514           2507.5816           2571.5682           6004.5853 
X181_HARV_ndvi_crop X197_HARV_ndvi_crop X213_HARV_ndvi_crop X229_HARV_ndvi_crop 
          8780.7595           8941.2789           8783.9587           8818.7488 
X245_HARV_ndvi_crop X261_HARV_ndvi_crop X277_HARV_ndvi_crop X293_HARV_ndvi_crop 
          8503.5736           7955.7723            332.0345            561.2631 
X309_HARV_ndvi_crop 
          5457.6267 

tower_ndvi <- raster::extract(ndvi_Z18, tower)
tower_ndvi

     X005_HARV_ndvi_crop X037_HARV_ndvi_crop X085_HARV_ndvi_crop X133_HARV_ndvi_crop
[1,]            3253.704            1888.675            2631.495            5887.102
     X181_HARV_ndvi_crop X197_HARV_ndvi_crop X213_HARV_ndvi_crop X229_HARV_ndvi_crop
[1,]            8793.096            8942.353            8867.162            8856.346
     X245_HARV_ndvi_crop X261_HARV_ndvi_crop X277_HARV_ndvi_crop X293_HARV_ndvi_crop
[1,]            8583.951            7804.457             331.571            563.1186
     X309_HARV_ndvi_crop
[1,]            5366.162

ndvi_df <- as.data.frame(tower_ndvi[1,])
ndvi_df

                    tower_ndvi[1, ]
X005_HARV_ndvi_crop       3253.7035
X037_HARV_ndvi_crop       1888.6754
X085_HARV_ndvi_crop       2631.4945
X133_HARV_ndvi_crop       5887.1022
X181_HARV_ndvi_crop       8793.0962
X197_HARV_ndvi_crop       8942.3530
X213_HARV_ndvi_crop       8867.1619
X229_HARV_ndvi_crop       8856.3461
X245_HARV_ndvi_crop       8583.9509
X261_HARV_ndvi_crop       7804.4574
X277_HARV_ndvi_crop        331.5710
X293_HARV_ndvi_crop        563.1186
X309_HARV_ndvi_crop       5366.1617

julian_day <- gsub("X|_HARV_ndvi_crop", "", row.names(ndvi_df))
julian_day

 [1] "005" "037" "085" "133" "181" "197" "213" "229" "245" "261" "277" "293" "309"

origin <- as.Date("2011-01-01")
ndvi_df$date <- origin + as.integer(julian_day)
names(ndvi_df) <- c("tower_ndvi", "date")
ndvi_df

                    tower_ndvi       date
X005_HARV_ndvi_crop  3253.7035 2011-01-06
X037_HARV_ndvi_crop  1888.6754 2011-02-07
X085_HARV_ndvi_crop  2631.4945 2011-03-27
X133_HARV_ndvi_crop  5887.1022 2011-05-14
X181_HARV_ndvi_crop  8793.0962 2011-07-01
X197_HARV_ndvi_crop  8942.3530 2011-07-17
X213_HARV_ndvi_crop  8867.1619 2011-08-02
X229_HARV_ndvi_crop  8856.3461 2011-08-18
X245_HARV_ndvi_crop  8583.9509 2011-09-03
X261_HARV_ndvi_crop  7804.4574 2011-09-19
X277_HARV_ndvi_crop   331.5710 2011-10-05
X293_HARV_ndvi_crop   563.1186 2011-10-21
X309_HARV_ndvi_crop  5366.1617 2011-11-06

ggplot() + 
  geom_point(data = ndvi_df, aes(x = date, y = tower_ndvi)) 

field_locations <- read.csv("./data/HARV_PlotLocations.csv")

class(field_locations)

[1] "data.frame"

head(field_locations)

   easting northing geodeticDa utmZone   plotID stateProvi    county domainName
1 731405.3  4713456      WGS84     18N HARV_015         MA Worcester  Northeast
2 731934.3  4713415      WGS84     18N HARV_033         MA Worcester  Northeast
3 731754.3  4713115      WGS84     18N HARV_034         MA Worcester  Northeast
4 731724.3  4713595      WGS84     18N HARV_035         MA Worcester  Northeast
5 732125.3  4713846      WGS84     18N HARV_036         MA Worcester  Northeast
6 731634.3  4713295      WGS84     18N HARV_037         MA Worcester  Northeast
  domainID siteID    plotType  subtype plotSize elevation  soilTypeOr plotdim_m
1      D01   HARV distributed basePlot     1600    331.64 Inceptisols        40
2      D01   HARV       tower basePlot     1600    341.62 Inceptisols        40
3      D01   HARV       tower basePlot     1600    347.61 Inceptisols        40
4      D01   HARV       tower basePlot     1600    334.34   Histosols        40
5      D01   HARV       tower basePlot     1600    352.93 Inceptisols        40
6      D01   HARV       tower basePlot     1600    342.43   Histosols        40

utm18crs <- st_crs(tower)
utm18crs

Coordinate Reference System:
  EPSG: 32618 
  proj4string: "+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs"

field_locations_sf <- st_as_sf(field_locations, coords = c("easting", "northing"),
                               crs=utm18crs)

field_locations_sf

Simple feature collection with 21 features and 14 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 731405.3 ymin: 4712845 xmax: 732275.3 ymax: 4713846
epsg (SRID):    32618
proj4string:    +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
First 10 features:
   geodeticDa utmZone   plotID stateProvi    county domainName domainID siteID
1       WGS84     18N HARV_015         MA Worcester  Northeast      D01   HARV
2       WGS84     18N HARV_033         MA Worcester  Northeast      D01   HARV
3       WGS84     18N HARV_034         MA Worcester  Northeast      D01   HARV
4       WGS84     18N HARV_035         MA Worcester  Northeast      D01   HARV
5       WGS84     18N HARV_036         MA Worcester  Northeast      D01   HARV
6       WGS84     18N HARV_037         MA Worcester  Northeast      D01   HARV
7       WGS84     18N HARV_038         MA Worcester  Northeast      D01   HARV
8       WGS84     18N HARV_039         MA Worcester  Northeast      D01   HARV
9       WGS84     18N HARV_040         MA Worcester  Northeast      D01   HARV
10      WGS84     18N HARV_041         MA Worcester  Northeast      D01   HARV
      plotType  subtype plotSize elevation  soilTypeOr plotdim_m
1  distributed basePlot     1600    331.64 Inceptisols        40
2        tower basePlot     1600    341.62 Inceptisols        40
3        tower basePlot     1600    347.61 Inceptisols        40
4        tower basePlot     1600    334.34   Histosols        40
5        tower basePlot     1600    352.93 Inceptisols        40
6        tower basePlot     1600    342.43   Histosols        40
7        tower basePlot     1600    343.11   Histosols        40
8        tower basePlot     1600    353.55 Inceptisols        40
9        tower basePlot     1600    338.11 Inceptisols        40
10       tower basePlot     1600    343.74   Histosols        40
                   geometry
1  POINT (731405.3 4713456)
2  POINT (731934.3 4713415)
3  POINT (731754.3 4713115)
4  POINT (731724.3 4713595)
5  POINT (732125.3 4713846)
6  POINT (731634.3 4713295)
7  POINT (731754.3 4713295)
8  POINT (731784.3 4712845)
9  POINT (731904.3 4713595)
10 POINT (731514.3 4713235)

st_write(field_locations_sf, "./data/FieldLocations_HARV.shp", 
         driver = "ESRI Shapefile", delete_dsn=TRUE)

Deleting source `./data/FieldLocations_HARV.shp' using driver `ESRI Shapefile'
Writing layer `FieldLocations_HARV' to data source `./data/FieldLocations_HARV.shp' using driver `ESRI Shapefile'
Writing 21 features with 14 fields and geometry type Point.